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Abstract 


Long-term, global offline (land-only) simulations with a dynamie vegetation phenology 
model are used to examine the control of hydroclimate over vegetation-related quantities. First, 
with a control simulation, the model is shown to capture successfully (though with some bias) 
key observed relationships between hydroclimate and the spatial and temporal variations of 
phonological expression. In subsequent simulations, the model shows that: (i) the global spatial 
variation of seasonal phenological maxima is controlled mostly by hydroclimate, irrespective of 
distributions in vegetation type, (ii) the occurrence of high interannual moisture-related 
phenological variability in grassland areas is determined by hydroclimate rather than by the 
specific properties of grassland, and (hi) hydroclimatic means and variability have a 
corresponding impact on the spatial and temporal distributions of gross primary productivity 
(GPP). 
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1, Introduction 


Recognition that the Earth’s energy and water cycles are intrinsically entwined is 
longstanding (e.g., Budyko 1971). The land surface energy and water balances both feature 
evapotranspiration as a dominant term, and the generation of rainfall (a key component of the 
water cycle) has a profound effect on the heat budget of the atmosphere. The inseparability of 
the energy and water cycles underlies their joint treatment in numerous analyses (e.g., Trenberth 
et al. 2011) and the formation of international research projects addressing their linkage, such as 
GEWEX (the Global Energy and Water Exchanges Project, part of the World Climate Research 
Programme, or WCRP). 

The Earth’s carbon cycle is in turn intrinsically entwined with the energy and water 
cycles. Vegetation health (and associated carbon uptake) is affected by water availability; 
deserts, for example, tend not to be carbon sinks. Conversely, carbon affects the water and 
energy cycles; the transpiration of water from vegetation and the associated cooling of the land 
surface are in large part controlled by the efficiency of the vegetation’s uptake of carbon dioxide 
(e.g.. Berry et al. 2010), and the build-up of vegetation through carbon uptake has a direct impact 
on land surface albedo - how much of the sun’s radiation is absorbed by the surface. Carbon 
dioxide is, of course, also a greenhouse gas. The basic connection between the surface fluxes of 
water, energy, and carbon is appropriately recognized in numerous studies (e.g., Eeuning et al. 
2004; Bowling et al. 2010), and it is a motivation for such international research projects as 
lEEAPS (the Integrated Eand Ecosystem Atmosphere Study, another component of WRCP). 

In this paper, we focus in particular on the carbon-water linkage at the land surface. A 
number of relevant studies in the literature have shared this focus. Using data collected at a 
number of flux tower sites in North America, Knapp and Smith (2001) provided a powerful. 
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geographically diverse analysis of the connections between the surface water and carbon cycles - 
specifically, of the controls of precipitation means and variability on aboveground net primary 
production (ANPP). Their results show that carbon uptake by the land surface is indeed strongly 
regulated by precipitation characteristics, with maximum uptake related strongly to precipitation 
amount and with the interannual variability of the uptake maximized in grassland areas, where 
both precipitation variability and vegetation cover are adequately high. Remotely sensed 
measurements of vegetation properties allow for an even more comprehensive and large-scale 
analysis of connections between carbon and climatic variables, including precipitation (e.g.. 

Fang et al. 2005, lichi et al. 2010, Jahan and Gan 2011). In a recent global analysis, Zeng et al. 
(2013) uncovered strong relationships between the interannual variations contained in a 
multidecadal NDVI dataset (normalized difference vegetation index, an indicator of green leaf 
area) and antecedent precipitation levels, particularly in temperate and tropical grasslands. 

A modeling framework is a natural venue for studying the connections between carbon 
and water. Wang and Eltahir (2000), using a simple coupled biosphere-atmosphere model, 
showed how the interaction between vegetation and precipitation can lead to multiple equilibria 
for vegetation state. Zeng et al. (1999) showed, again with a simple coupled model, how 
vegetation-climate interactions may affect the nature of precipitation variability in the Sahel. 
Puma et al. (2013) used a modeling framework to compare the impacts of meteorological 
variability and phenological variability on the simulation of surface moisture and carbon fluxes. 
Complex and relatively complete models of vegetation behavior, models that indeed tie together 
explicitly the interactions between carbon, energy, and water fluxes at the land surface and 
accordingly allow the prediction of vegetation state, are arguably the new state-of-the-art in 
numerical climate modeling. Sellers et al. (1997) pointed to the explicit treatment of carbon as a 
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logical step in the evolution of land surfaee treatments in Earth system models; dynamic 
vegetation models (DVMs) following this evolutionary path are already being used at major 
climate modeling centers (e.g., Lawrenee et al. 2010, Krinner et al. 2005, Boussetta et al. 2012, 
Dunne et al. 2013). 

An advantage of using a modeling framework for earbon-water studies is the potential for 
doing unique analyses that isolate and illustrate the mechanisms that control the transfers of 
water and earbon aeross the land surfaee. Carefully formulated modifieations of a physical 
process treatment or of a variable that forces it can be imposed, and the resulting impacts on 
surface fluxes can be quantified and analyzed, thereby elueidating the role of the proeess 
examined. A second important advantage of such models is their ability to provide data fields 
that are unattainable with in situ measurement networks or even satellite -based sensors. Gross 
primary productivity (GPP), for example, can only be measured direetly at a limited number of 
flux tower sites. A DVM, however, if driven with observations-based meteorologieal foreing, 
can potentially produce estimates of GPP at high spatial and temporal resolution across the 
globe. Sueh estimates would be biased relative to nature, of eourse, due to deficieneies in model 
formulation and foreing data; still, if care is given to their interpretation, the estimates do have 
seientific value. 

Both of these advantages eome into play in the present paper, in whieh we use the 
dynamic phenology eomponent of an established DVM together with the water and energy 
balance framework of a hydrology-focused LSM to eharaeterize, on a global seale, the controls 
of precipitation means and variability on GPP - both on its spatial distribution and on its 
temporal variability across the globe. The modeling system used (described in Seetion 2) is 
indeed found to be effective in capturing the key hydroclimatic eontrols on phenology that 
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operate in nature (as demonstrated in Seetion 3). The simulated GPP distributions from the thus- 
validated system are analyzed jointly with global preeipitation data in Seetion 4. The model 
experiments provide new insights into the relative impaets of precipitation means, precipitation 
variability, and vegetation type in determining GPP distributions. 


2, Dynamic Phenology Model 

The dynamic phenology model used in this study is in essence a merger of the carbon 
physics of the NCAR/DOE CLM4 dynamic vegetation model (Oleson et al. 2010) with the 
energy and water balance formulations of the NASA Global Modeling and Assimilation Office 
(GMAO) Catchment land surface model (LSM) (Koster et al. 2000). We provide here a brief 
description of these two components and the technique used to merge them into a new model, 
hereafter referred to as the Catchment-CN LSM (i.e., the Catchment LSM with carbon and 
nitrogen physics). 

The NCAR/DOL Community Land Model, version 4 (CLM4), represents prognostic 
coupled energy, water, carbon, and nitrogen cycles in a framework that permits global-scale as 
well as regional and site-level simulation. The global-scale parameterization used here includes 
specification of sub-grid heterogeneity in plant functional type (PLT) distributions, with multiple 
PLTs assigned fractional area coverage within each grid cell, where they compete with one 
another for available soil moisture and mineral nitrogen resources. In this prescribed 
biogeography mode the fractional areas occupied by individual PLTs do not change, but 
vegetation growth, soil heterotrophic activity, carbon stocks, and other ecosystem states (such as 
leaf area index) do vary prognostically (Thornton et al. 2009). 
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The GMAO Catchment land surface model is a state-of-the-art surface energy and water 
budget model designed for use with global Earth system models. As with most other LSMs, the 
Catchment LSM employs complex treatments of land surface flux generation, tying the 
efficiency of evaporation and runoff generation to the moisture and temperature states of the land 
surface, and it includes parameterizations of vegetation impacts on transpiration, canopy 
interception, albedo, and surface roughness. Relatively unique to the Catchment LSM is its 
treatment of the subgrid variability of soil moisture and temperature, which is explicitly tied to a 
description of the topographic variability in the region modeled - in the Catchment LSM, valley 
bottoms within a given grid element are explicitly modeled as being wetter, and the hilltops are 
explicitly modeled as being drier. Runoff and evaporation are calculated independently in the 
different hydrological regimes, using regime-specific physics. 

In essence, in merging the two models, we retain the Catchment LSM’s energy and water 
balance calculation framework while using the NCAR/DOE CLM4 carbon balance calculations. 
The approach is illustrated in Ligure 1. In the original Catchment LSM (Ligure la), the model 
uses forcing from the atmosphere along with prescribed vegetation phenology (LAI and 
greenness fraction) and the current values of LSM temperature and moisture prognostic variables 
to compute the canopy conductance, the parameter describing the ease with which the plants 
transpire water. The canopy conductance, computed separately for each hydrological regime, is 
then used in each regime’s energy balance and water balance calculations, which in turn provide 
the fluxes of heat and moisture to the atmosphere. 

Ligure lb shows the approach used by the merged system, the Catchment-CN LSM. The 
atmospheric inputs are now fed first into the components of the NCAR/DOE model that update 
the carbon states and compute, as a matter of course, canopy conductances that reflect an explicit 
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treatment of photosynthesis physics. These canopy conductances, along with the leaf area 
indices diagnosed from the new carbon prognostic variables, are fed into the energy and water 
balance calculations of the original Catchment LSM. The output fluxes with the merged system 
include a net carbon flux. 

The merger of the two models allows the Catchment-CN LSM to follow 1 9 distinct 
vegetation types, a significant increase from the six independent types followed with the original 
Catchment LSM. Furthermore, the unique character of the original Catchment LSM allows for 
the independent monitoring of carbon variables in the different topographically-defined 
hydrological regimes. Figure 2 describes our methodology. Each land surface element is 
subdivided into three static carbon zones defined by topography, through analysis of the 
distribution of the compound topographic index (Moore et ah, 1993). The first zone, covering a 
fixed 10% of the area, represents the valley bottoms; this zone tends to be generally wet. The 
second and third zones represent the lower (drier) hillslopes and upper (even drier) hillslopes, 
respectively. Through areal weighting, soil moisture and temperature information from the 
dynamically-varying hydrological zones are combined for use by the carbon physics in the fixed 
vegetation zones, as indicated in the figure. Separate sets of carbon prognostic variables are 
followed in each vegetation zone, and thus each zone generates a different manifestation of 
phenology. When examining the model results, we find that green vegetation indeed tends to be 
densest in the valley bottoms. 

Some additional modifications to the NCAR/DOE vegetation model were needed to 
optimize its performance in the GMAO system. To prevent some occasional singular behavior - 
namely, the catastrophic shutdown of vegetation during cold spells and a resulting overgrowth of 
the vegetation during the subsequent year - we replaced a particular set of vegetation types (crop 
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and temperate shrubs/grass) that feature a strong response to temperature stress by a mix of two 
different types: one that is seasonally deeiduous and one that is not. Neither of the replaeement 
types employ the temperature stress shutdown, though both respond to moisture stress; the 
proportion of the mix applied is defined by latitude, and the replaeement is indeed limited to the 
latitude band 32 °- 42 ° in both hemispheres. Outside of this latitudinal band, we limit the number 
of eoexisting PFTs in eaeh statie earbon zone to two. Also, we modified the NCAR/DOE 
vegetation physios to allow half of the new earbon assimilated by deeiduous types to be 
displayed during the ourrent year rather than in the following year, whioh brings oertain 
measures of our interannually-varying phenology more in line with observations. Finally, 
whereas the NCAR/DOE vegetation model uses the previous year’s annual mean temperature to 
determine oertain onset triggers, we use a olimatologioal mean temperature. 

In our main (“oontrof ’) applioation of the model, the presoribed distributions of 
vegetation type follow those used by the default 0.5°x0.5° version of CLM4 (Oleson et al. 2010). 
Vegetation phenology and earbon states, however, evolve freely. The model is run globally 
offline (i.e., disoonneoted from an atmospherio model) on high-resolution eatehments (roughly 
20-30 km in size) over the period 1948-2008, using the observations-based meteorologieal 
foreing of Sheffield et al. (2006); the simulation loops over this period more than 30 times to 
ensure spin-up and equilibration of the earbon storage reservoirs. The output data examined 
(phenologieal variables, earbon fluxes, ete.) are aggregated to 2°x2.5° for proeessing. 


3. Evaluation Against Observations 

To test the realism of the model’s connections between hydroclimate and vegetation 
variables, we focus on two distinct aspects of global phenologieal expression: the global spatial 
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pattern of long-term phenologieal means and the interannual variability of phenology at a given 
loeation. These are diseussed in turn following a brief deseription of the observations. 


a. Observations used 

We examine satellite-based produets of NDVI (normalized differenee vegetation index) 
and FPAR (fraetion of absorbed photosynthetieally aetive radiation), both of whieh inerease with 
green vegetation eover. The NDVI data is a subset of the latest version of the Global Inventory 
and Mapping Studies, or GIMMS, data (Tueker et al. 2005). The data’s native resolution is 
semiweekly at 8 km, and the data span the period July 1981 -present. For our analyses we 
aggregate these data to a 2.5°x2.5° degree, monthly resolution for the period 1982-2010. The 
data are derived from the Advaneed Very High Resolution Radiometer (AVHRR) instrument 
with known limitations eompared to the more advaneed MODIS instrument (Kaufman et al. 
1998). However, the longer temporal eoverage of GIMMS relative to MODIS (29 versus 1 1 
years) and the good eorrespondenee between their measurements (Tueker et al. 2005, Beek et al. 
2011) makes it well suited to the analysis presented here. 

The FPAR data are derived direetly from the NDVI data using the method of Los et al. 
(2000). The method eombines the NDVI-based FPAR estimation teehnique of Sellers et al. 
(1996) with that of Choudhury (1987) and Goward and Huemmrieh (1992); the eombination 
provides estimates that are well behaved relative to available in-situ observations. The 
relationship between NDVI and FPAR underlying this eombined approach is monotonic but 
nonlinear. Note that it is also somewhat vegetation dependent, so that the conversion of global 
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NDVI data to global FPAR data requires a global field of vegetation types. Thirteen years of 
FPAR data are available, spanning the period 1997-2009. 

As will be seen below, the sensitivities of the NDVI and FPAR data to hydroclimatie 
variation are similar in many ways. Both are worth illustrating here. The NDVI values are 
eonstructed direetly from speetral reflectanee measurements and thus represent a raw form of the 
observations. While the construetion of the FPAR values requires some additional assumptions 
regarding vegetation behavior, FPAR has the distinet advantage of representing a physieally 
meaningful phenologieal variable, one that ean be compared directly to output from the 
Catchment-CN model. 

The global precipitation data used here consist of monthly precipitation totals for 1979- 
present at 2.5°x2.5° degree resolution, as produced by the Global Precipitation Climatology 
Project as part of their Version 2 Satellite-Gauge dataset (Adler et al. 2003; see also 
ftp://precip.gsfc.nasa.gOv/pub/gpcp-v2.2/docW2.2 doc) . Satellite -based data contributing to the 
product, in varying capacities and over various periods and regions, include Special Sensor 
Microwave/Imager (SSM/I) passive microwave estimates, Television-Infrared Observation 
Satellite (TIROS) Operational Vertical Sounder (TOVS) estimates, and the Adjusted 
Geostationary Operational Environmental Satellite (GOES) Precipitation Index (Adler et al. 
1994). A wealth of surface rain gauges is used to adjust the multi-sensor precipitation estimates 
over land. Hall et al. (2006) provide background on the accuracy of the GPCP product; of note is 
the higher uncertainty of the product over mountains, deserts, high latitudes, and undeveloped 
areas due in large part to a lower density of rain gauges. 
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b. Impact of hydrological variations on the mean spatial distributions of phonological variables 

To deal with the faet that NDVI shows significant seasonal variability, with different 
regions having different peak months for the index, we examine a quantity we will eall NDVImax- 
We eompute, at eaeh 2.5°x2.5° grid cell, the average seasonal cyele of NDVI from the GIMMS 
data and then identify the month for whieh the average NDVI is highest. NDVImax is set to the 
average value for the 3 -month period centered on this peak month. (Note that under this 
definition, the values for NDVImax in adjaeent grid cells may be taken from different 3-month 
periods.) Figure 3a shows the global distribution of NDVImax as derived from the spatially 
aggregated GIMMS data. The distribution mirrors known vegetation distributions, with large 
values in tropical, deciduous, and boreal forests, intermediate values in grassland and shrubland 
areas, and small values in the deserts. 

Figure 4a shows how the spatial distribution of NDVImax in Figure 3a eorrelates with 
various meteorologieal quantities. The first four bars of eaeh panel show the square of the 
spatial correlation (r , aeross land surface grid cells) of NDVImax with, respectively, annual mean 
preeipitation, the standard deviation of annual preeipitation, annual mean air temperature, and 
annual mean net radiation. (The base- 10 logarithms of the preeipitation quantities are in faet 
used here. Temperature and net radiation information are derived from the full period of the 
Sheffield et al. (2006) dataset. The outgoing longwave eomponent of the net radiation is 
estimated using the surfaee air temperature in that dataset.) The salient result from the figure is 
the dominanee of the two preeipitation quantities in determining the spatial structure of NDVI. 
Multiple regression of NDVImax on the mean and variability of preeipitation produees an of 
about 0.55 (fifth bar), and adding the temperature and net radiation information to the multiple 
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regression does not signifieantly inerease the (sixth bar). These results underlie the 
importanee of hydroelimate in determining the spatial distribution of phenologieal maxima. 

Figure 5a shows more direetly how NDVImax is related to preeipitation means and 
variability. Eaeh dot in Figure 5a eorresponds to a 2.5°x2.5° land grid cell. The size and color 
of the dot is determined by the local value of NDVImax, as indicated by the legend. The dot’s 
abscissa is determined by the mean annual precipitation at that grid cell, and the dot’s ordinate 
refers to the interannual variability of precipitation there. (Note the logarithmic scales.) The 
precipitation and NDVI quantities are computed over consistent time periods; for example, if a 
grid cell’s peak NDVI, as computed from the GIMMS data for 1982-2010, occurs in July, then 
precipitation means and variances are computed from nineteen September- August yearly totals 
starting with the total for the period September 1981 -August 1982. 

Two features of the scatter plot stand out. The first reflects an expected result: a 
minimum average precipitation must be achieved to attain moderately high NDVImax levels. The 
plot shows this minimum value to be roughly 1 mm/day; the dots to the left of this threshold 
(which include, of course, all desert points) show low values of NDVImax- The second, and more 
intriguing, feature of the scatter plot is the tendency for NDVImax to decrease as the standard 
deviation of precipitation increases. This feature is illustrated more clearly in Figure 5b, which 
shows a binned version of the scatter plot data; to generate this plot, an array of boxes is overlain 
on Figure 5a, and the NDVImax values for the points within each box are averaged. For a given 
value of the mean precipitation, especially for values above 1 mm/day, NDVImax clearly tends to 
decrease with increasing Gp. This presumably reflects the reduced ability of vegetation to 
flourish when the year-to-year supply of water is less stable. 
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We also examine in this context the analogous variable FPARmax, the average value of a 
grid cell’s FPAR for the 3-month period centered on the peak FPAR month, as determined from 
the local climatological cycle. Figure 3b shows the distribution of FPARmax as computed from 
the GIMMS data. As might be expected, given that FPAR in GIMMS is derived from NDVI, the 
spatial distributions in Figures 3 a and 3b are very similar, as are the spatial correlations with the 
meteorological forcing variables (Figure 4b). Figures 5c and 5d show the precipitation-based 
scatter plots for the FPARmax values. Average water supply (mean precipitation) and water 
supply stability (a?) are seen to impose dual control over FPARmax as well; the sensitivity of 
FPAR to hydroclimate is indeed very similar to that of NDVI. 

How well does the Catchment-CN model perform? The model produces diagnostics for 
both the incident and absorbed photosynthetically active radiation; we take the ratio of these 
quantities to produce the model’s FPAR values. Figure 3c shows the global distribution of 
simulated FPAR in the peak 3 -month period; note that for a given location, this peak period may 
differ from that for the observations. Two features of the simulated FPAR distribution stand out. 
First, the simulated spatial patterns agree well with the observed patterns in Figure 3b. Second, 
there are, nonetheless, apparent biases in the simulated FPAR values, with the highest simulated 
values being too large and the lowest being too small. Such biases presumably reflect 
deficiencies in the model, though they may also stem partially from limitations in the forcing 
data or in the observational FPAR values themselves. The biases must be kept in mind 
throughout our analysis. 

The square of the spatial correlation of simulated FPAR with meteorological forcing 
variables (Figure 4c) agrees quite well with the corresponding values found for observed FPAR 
(Figure 4b). The simulated r values with the temperature and net radiation variables are slightly 
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higher, but these values are still quite low. Figure 4 shows that, in strong agreement with the 
observations, variations in hydroelimate explain most of the FPAR variability seen in the model. 

The agreement in spatial pattern with a presenee of bias also manifests itself in the 
preeipitation-based seatter plot in Figure 5e and the eorresponding binned version of the plot in 
Figure 5f. In agreement with the observations, the model elearly shows an inerease in FPAR 
with inereasing preeipitation and with deereasing preeipitation variability. Overall, the model, 
though biased, does appear to simulate realistie eontrols of hydroelimatie variation over 
phenologieal means. 


e. Impaet of hydrologieal variations on the interannual variability of phenologieal variables 

As a seeond and somewhat independent test of the ability of the Catehment-CN model to 
eapture observed links between earbon and water variables, we examine the interannual 
variability of vegetation phenology. Rather than examining the total varianee of a variable sueh 
as summertime NDVI, we foeus instead on a modified quantity, one that eaptures the earbon- 
water eonneetion: 

Var(NDVI)* = Var(NDVI) Corr^(NDVI,P) , (1) 

where Var(NDVI) is the interannual varianee of 3-month NDVI averages (again eentered on the 
peak NDVI month, based on the elimatologieal seasonal eyele), Corr (NDVI,P) is the eorrelation 
between these individual NDVI averages and the eorresponding yearly preeipitation totals (with 
the end of the preeipitation averaging period eorresponding to the end of the 3-month NDVI 
averaging period), and Var(NDVI)* is interpreted as the portion of the NDVI varianee assoeiated 
with variations in moisture availability. That is, we are employing here the standard 
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2 

interpretation of Corr (NDVI,P) as the fraetion of the varianee of NDVI “explained” by 
variations in P. Equation (1) allows us to isolate this part of NDVI variability from that 
associated with other sources, such as variations in radiation or nutrients as well as interference 
from clouds, water vapor, and aerosols (Los et al. 2000). 

A few notes are required regarding the estimation of Var(NDVI)*. First, by using the 
annual totals for precipitation, we are assuming that a given year’s precipitation represents the 
water available that year for growth. Of course, other averaging periods for the precipitation 
could have been employed (e.g., Zeng et al. 2013). The patterns in Corr (NDVI,P) obtained with 
these other averaging periods, however, turn out to be the same, to first order; correlation maps 
generated using 6-month or 9-month precipitation averages (not shown) are very similar to those 
generated with the annual precipitation. Note that using the annual precipitation rather than the 
contemporaneous 3-month precipitation has an important advantage: it reflects the fact that 
antecedent precipitation can provide water to vegetation growth through storage in ground 
reservoirs and snowpack (Milly, 1994). 

Second, the observations are known to be subject to significant contamination from 
clouds in high latitudes and from pollution in Southeast Asia (Fensholt and Proud, 2012), the 
upshot being that small and artifactual negative correlations between NDVI and precipitation are 
often seen in these regions. These negative correlations are problematic for our analysis. We 
zero them out before computing Corr (NDVI,P), making the explicit assumption that any such 
negative correlations represent noise. Note that even on the off chance that the negative 
correlations are real, they would not represent the physical relationship we are after in this paper, 
namely, the ability of water limitations to limit vegetation growth. 
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Figure 6a shows the distribution of Var(NDVI)*, as eomputed with (1). The patterns are 
quite interesting: the regions for whieh moisture-related NDVI variability is high tend to 
coincide with the Earth’s grassland regimes - in the Great Plains of the U.S., the Nordeste region 
of Brazil, the African Sahel, the Asian steppes, and eastern and northern Australia (see Figure 7). 
The Var(NDVI)* patterns do miss grassland areas in India and China, but as shown in Figure 7, 
these areas are subject to extensive irrigation (Siebert et al. 2005), a supply of water not 
accounted for in the Corr (NDVI,P) diagnostic. Figure 6 demonstrates that, aside from such 
irrigated areas, the locations of the Earth’s grassland areas can be identified reasonably well from 
the joint analysis of NDVI and precipitation data. The same patterns, and thus the same 
connections to grassland regimes, are seen for Var(EPAR)*, the portion of the interannual 
variance in 3-month EPAR averages related to moisture variations. 

The results obtained with the dynamic phenology model are remarkably similar. A 
comparison of Eigures 6b and 6c shows that the model captures very well the observed spatial 
pattern Var(EPAR)*, though again with a bias, as indicated by the different scaling factors used 
for the plotting. Overall, the model successfully captures the role of hydroclimate in 
determining the spatial distribution of interannual variability in phenology. 


4, Model Experiments 

Having demonstrated the Catchment-CN model’s ability to capture the basic 
hydroclimatic controls on phenology seen in the observations, we now use model experiments to 
address key questions regarding the connections between hydroclimate and vegetation. 
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a. Influence of vegetation type on phonological variability 


Clouding the interpretation of the Catchment-CN model’s performance relative to 
observations in Section 3 above is the possibility that its use of prescribed vegetation types is 
somehow guaranteeing correct model behavior. Given, for example, that the observed 
distribution of Var(FPAR)* in Figure 6b captures well the locations of the world’s grasslands 
(Figure 7), we must consider the possibility that high values of Var(FPAR)* are encouraged by 
the unique properties of grassland and discouraged by the properties of forests and shrubs, so that 
by imposing the observed vegetation distributions in the model, we artificially guarantee high 
simulated values of Var(FPAR)* in the correct areas (Figure 6c). The more intriguing possibility 
to consider, however, is that a specific hydroclimatic regime is responsible for high Var(FPAR)* 
values, a regime for which only grasslands happen to survive. With this second possibility, the 
vegetation type does not cause the Var(FPAR)* value; rather, the vegetation type and the 
Var(FPAR)* value are together controlled by something else, namely, the local moments of 
precipitation. 

To examine this issue, we performed a repeat of the simulation described above, but with 
a twist: grassland vegetation was imposed on all land surfaces, and no other vegetation types 
were allowed to exist. Thus, in this experiment, vegetation type could not affect in any way the 
simulated spatial and temporal distributions of FPAR. Note that in this experiment, grassland is 
placed even in the driest deserts and in the wettest tropical areas; if the local climate is not 
conducive to grassland’s survival, the grass is accordingly allowed to die out. 

Figure 8c shows the spatial distribution of Var(FPAR)* for the all-grassland simulation. 
The plot captures, to first order, the features seen in the original model plot, supporting the 
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second possibility noted above. That is, the presence of grassland does not lead to high 
Var(FPAR)* values; the high values are instead indicative of a hydroclimatic regime that also 
happens to support grassland best. Similarly, the all-grassland simulation shows a relationship 
between FPAR maxima, mean precipitation, and precipitation variability (Figure 8a) that agrees 
to first order with that seen in the original model simulation (Figure 5f). The fact that FPAR 
tends to be highest in very wet conditions, for example, is not simply the result of the presence of 
dense forests in wet areas; the wet conditions themselves encourage the high FPAR values, and 
wet areas also tend to be where dense forests tend to flourish. 

We repeated the simulation still again, this time after prescribing a deciduous forest 
vegetation type everywhere. The results, shown in Figures 8b and 8d, are essentially the same. 
Hydroclimatic variability, more than vegetation type, appears to dominate phenological 
variability - in the model and, we can infer, in nature. 


b. Hydroclimate and the global carbon cycle 

As noted in the introduction, a unique advantage of a model that can simulate phenology 
is its ability to provide information on additional, difficult to measure quantities. While carbon 
fluxes such as gross primary productivity (GPP), net primary productivity (NPP), and net 
ecosystem exchange (NEE) have been measured at various tower sites (Baldocchi 2008), directly 
observed global distributions of land-atmosphere carbon exchange are nonexistent. Model 
simulations, however, can readily provide these fields, and many examples of such simulated 
distributions already appear in the literature (e.g., Eriedlingstein 2006). (We note that other 
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approaches for inferring global fields, sueh as maehine learning algorithms that upseale from the 
site measurements, are also available [Jung et al. 2011].) 

Here we provide model-based estimates of the eonneetion between earbon exehange and 
hydroelimatie variability, foeusing mainly on GPP. We first provide in Figure 9a this partieular 
model’s vision of the global distribution of GPP. Beeause GPP is a flux rather than a 
manifestation of vegetation state, we present it in terms of annual averages rather than for a 3- 
month maximum period. The distributions have the expeeted maxima in the densely forested 
tropies, with swaths of high values in the boreal forests of the north. Figures 9b and 9e show the 
eorresponding GPP fields from the simulations preseribing grassland and deeiduous tree types, 
respeetively. The three panels show some differenees but are, to first order, very similar, 
indieating that vegetation type alone is not the main souree of spatial variations in GPP; both 
GPP and vegetation distributions are apparently eontrolled in tandem by something else. 

Naturally, that “something else” is water availability. Figure 10 shows, in analogy to 
Figure 4, the square of the spatial eorrelation between GPP and various meteorologieal foreing 
variables. For all three simulations (eontrol, “all grass”, and “all trees”), preeipitation mean and 
variability have the dominant impaet on GPP, with an r of about 0.55 for the multiple regression 
of GPP on logioP and logioop. Adding in the annual temperature and net radiation information 
inereases the r to about 0.65. The faet that the r values do not inerease by mueh for the 
uniform vegetation experiments suggests onee again that variations in vegetation type do not by 
themselves eontribute signifieantly to spatial variations in GPP; the remaining unexplained 
varianee in Figure 10a presumably results from spatial variability in, for example, the seasonal 
eyeles and shorter-term temporal strueture of the foreing quantities. 
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Figure 1 1 shows how precipitation means and variability control the spatial distribution 
of GPP using scatter plots analogous to those shown in Figure 5. As with FPAR, GPP tends to 
increase with increasing moisture availability (x-axis) and decreasing interannual variability (y- 
axis), regardless of which vegetation types are assigned at the surface. 

In contrast to Figure 5, Figure 1 1 uses a nonlinear scale for the shading, a scale that 
shows the dominance of precipitation means over precipitation variability in determining GPP. 
The impact of precipitation variability on GPP, however, is nevertheless significant. This is 
demonstrated with a supplemental model simulation (“ClimP”) in which we prescribed standard, 
spatially varying vegetation types (as in the control simulation) but a modified precipitation 
forcing: at each grid cell in ClimP, we scaled the precipitation forcing in each month of each year 
so that the seasonal cycle of monthly totals for the year matched the long-term (climatological) 
seasonal cycle. Thus, in ClimP, we artificially removed the monthly-scale year-to-year temporal 
variability in the precipitation forcing - at each grid cell, the mean precipitation applied was 
identical to that used in the control simulation, whereas the interannual variability of monthly 
precipitation was, by construction, set to zero. 

Figure 12 shows the difference between the mean annual GPP produced in ClimP and 
that in the control simulation. Regions with large positive differences appear in the southeast 
U.S., along the eastern coasts of South America and Australia, in the Indian subcontinent, in 
northeastern China, and in various other regions of South America and Africa. Negative 
differences do not appear anywhere. In effect. Figure 12 illustrates where GPP in the real world 
would be larger if the year-to-year precipitation supply were more dependable - i.e., where the 
interannual variability of precipitation holds down the land surface’s carbon uptake. Note, 
however, that human activities can mitigate the effects of this variability. India, southeast Asia, 
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and northeastern China in partieular are known to undergo extensive irrigation (Figure 8). 
Beeause irrigation is effeetively a means of providing a more dependable water supply, these 
partieular areas may, in the real world, be eapturing the larger GPP rates. 

With Figure 13, we foeus on the interannual variability of GPP at eaeh grid eell rather 
than on the spatial distribution of its mean. Figure 13a shows the varianee of annual GPP. 

Figure 13b shows the spatial distribution of Corr (GPP,P), where P is the annual preeipitation; 
that is. Figure 13b shows the fraetion of the total GPP variance that is associated with, or can be 
“explained by”, variations in annual water supply. The fractions are reasonably large across the 
globe, even in some areas considered to be not strongly water-stressed, such as the southeastern 
United States. In contrast, the fields of Corr^(GPP,T) and Corr^(GPP,Rnet), where T is the 
yearly-averaged temperature and R„et is the yearly-averaged net radiation, show significantly 
lower values (Figures 13c and 13d). While interannual temperature variations do have some 
impact on high latitude GPP variations (perhaps through their effects on snowcover duration), 
they have little impact anywhere else. Interannual net radiation variations appear to contribute 
more, especially in Africa; it is quite possible, however, that these particular “contributions” are 
not real and instead simply reflect known existing correlations between precipitation and net 
radiation there (not shown). 

Together, annual precipitation, temperature, and net radiation do not explain all of the 
simulated GPP variability. As before, presumably a significant part of the variability stems from 
year-to-year variations in (for example) the sub-annual timing of the precipitation and associated 
variations in infiltration and runoff. 
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Figure 14 shows one final interesting result regarding the interannual variability of GPP. 
The shading shows Var(GPP) for a 3 -month averaging period (centered, at each grid cell, around 
the month of maximum GPP). Overlain on the plot are black dots indicating where Var(FPAR) 
for 3 -month averages (centered around the local monthly maximum for FPAR) exceeds a value 
of 0.003, an arbitrary threshold chosen for plotting convenience. The figure shows that Var(GPP) 
and Var(FPAR) tend not to be maximized in the same regions; Var(GPP) maxima tend to lie on 
the wetter sides of the Var(FPAR) maxima. The same basic result (not shown) is found for 
comparisons of the water-limited portions of the variances (i.e., Var(GPP)* versus Var(FPAR)*), 
and it is also found (not shown) for the all-grassland and all-deciduous-trees simulations, 
suggesting that variations in vegetation type are not responsible for such spatial offsets in the 
maxima. The spatial offsets are instead induced by the carbon physics built into the modeling 
system. Assuming these physical treatments are accurate, then similar offsets would apply to the 
real world’s distributions of Var(GPP) and Var(FPAR). In other words, given estimates of 
Var(FPAR) attained, for example, through the processing of the GIMMS data, knowledge of the 
offsets could potentially help in the construction of an estimated spatial field of Var(GPP). 


4. Summary and Discussion 

Using the Catchment-CN model, a merger of the dynamic phenology components of the 
CLM4 dynamic vegetation model with the water and energy budget framework of the GMAO 
Catchment LSM, we examine the connections across the globe between hydroclimate and 
vegetation variables. Justification for the use of this model in such a study is provided by its 
demonstrated ability to reproduce observed connections between FPAR and precipitation 
moments (Section 3), namely, the increase in FPAR with increasing mean precipitation and 
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decreasing precipitation variability and the proper geographical placement of spatial maxima in 
the global field of moisture-related FPAR variance. 

Our model results can be summarized as follows. First, based on our supplemental 
simulations with globally uniform vegetation type, we find that the aforementioned relationships 
between FPAR and precipitation moments are largely independent of vegetation type; the fact 
that trees grow in wet regimes, grass grows in drier regimes, and shrubs grow in even drier 
regimes has only a second-order impact on the spatial distribution of FPAR and its interannual 
variability at each location. Instead, hydroclimatic moments appear to be the dominant 
determinants of both vegetation type and phenological expression, as represented by FPAR. Our 
second basic result is that hydroclimatic moments provide a similarly dominant control over the 
spatial and temporal variability of gross primary productivity (GPP), again with only a second- 
order contribution from vegetation type. 

Such a global scale description of GPP connections to hydroclimate is achievable with a 
DVM but is not possible with observations, which are much more spatially and temporally 
limited. Knapp and Smith (2001) used observations collected across eleven tower sites to show 
that aboveground net primary production (ANPP) tends to increase with increasing annual 
precipitation, and our global scale results (for GPP, a related variable) are consistent with this. 
We do see some inconsistencies, however, with their study. For example, Knapp and Smith 
(2001) find that ANPP has its maximum interannual variability in grassland areas. We find that 
while the interannual variability of FPAR is maximized in grassland areas, the maxima for GPP 
variability tend to be spatially offset from these FPAR variance maxima (Figure 14), slightly 
toward the wetter (forested) side. The offset is minor, however, and the apparent inconsistency, 
while certainly a possible result of model deficiencies, may also relate to the limited number of 
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tower sites they examined. More importantly, Knapp and Smith (2001) find that “interannual 
variability in ANPP [is] not related to variability in preeipitation”. Results from our control 
simulation (not shown) indicate that the square of the spatial correlation coefficient between 
Var(GPP) and Var(P) across land points is of the order of 30%, which disagrees with their 
conclusion; indeed, when we limit the calculation to values at the grid cells containing the LTER 
sites they studied, the square of correlation coefficient increases. We also find a reasonably 
strong relationship between the time series of GPP and precipitation (Figure 13b) at individual 
locations. 

While interpretations of DVM -based results must be tempered by knowledge of model 
biases and limitations, DVM experiments, if properly interpreted, open the door to a wealth of 
potential studies of the global carbon cycle and its interactions with the global water and energy 
cycles. This paper provides one such study. Another example of note is provided by Guan et al. 
(2012), who show with DVM simulations over Africa that the statistical character of 
precipitation forcing (e.g., rainfall intensity) manifests itself in the GPP produced. The 
advantages of using DVMs - their provision of comprehensive (and often unmeasurable) data 
and their ability to be modified at will to allow the examination of the impacts of individual 
physical processes - stand them in good stead for future carbon analyses. Our understanding of 
global carbon-water-energy connections should continue to increase as researchers continue to 
use ever-improving versions of these tools. 
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Figure 1. a. Schematic of flux computations in the original Catchment LSM. b. Schematic of 
flux computations in the merged model, Catchment-CN. 

Figure 2. Schematic of independent vegetation (carbon) treatments in topographically-defined 
vegetation zones. Three static vegetation zones are defined, with independent carbon 
prognostic variables in each. Wi, W2, and W3 are soil moisture states in the three 
dynamically-varying hydrological zones (with time-varying areas ARl, AR2, and AR3); 
weighted averages of these states (e.g, Wvi, as shown in the figure) and corresponding 
weighed temperature states are passed down to the carbon physics calculations for the 
different vegetation zones. 

Figure 3. a. Distribution of average NDVI in peak NDVI season (the month for which the mean 
seasonal cycle of NDVI is maximized along with the preceding and following months), 
from GIMMS observations. The peak season varies with grid cell; see text for details, b. 
Same, but for average FPAR in peak FPAR season, from GIMMS observations, c. Same, 
but for average FPAR in peak FPAR season, from model simulation. 

Figure 4. (a) Square of the spatial correlation coefficient (over land grid cells) between observed 
NDVI and land surface forcing variables: the logarithm of the mean annual precipitation 
(P), the logarithm of the standard deviation of annual precipitation (op), mean annual air 
temperature (T), and mean annual net radiation (Rnet)- The final two bars show the square 
of the correlation coefficient obtained from the multiple regression of NDVI against, 
respectively, (i) precipitation mean and standard deviation, and (ii) all four quantities, (b) 
Same, but for observed FPAR. (c) Same, but for modeled FPAR. 
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Figure 5. (a) Average GIMMS NDVI in peak NDVI season as a funetion of the mean 

preeipitation (x-coordinate) and the standard deviation of annual preeipitation (y- 
coordinate). Eaeh dot represents a single land grid eell. (b) Same as (a), but with the 
individual values in the seatter plot averaged over bins. At least 5 dots must lie within a 
bin for the binned value to be plotted, (e) and (d): Same as (a) and (b), but for GIMSS 
FPAR data, (e) and (f): Same as (a) and (b) , but for model-simulated FPAR data. 

Figure 6. (a) Produet of the interannual varianee of GIMMS NDVI data averaged over the 
maximum NDVI season and the square of the eorrelation between NDVI and annual 
preeipitation, multiplied by 1000. (b) Same, but for GIMMS FPAR data, (e) Same, but 
for model-generated FPAR data, and with the sealing faetor ehanged to 500. 

Figure 7. Map of grassland loeations, as derived from the distributions used in the Seeond Phase 
of the Global Soil Wetness Projeet (Dirmeyer et al. 2006). The dots overlain on the plot 
indieate regions for whieh irrigation is extensive (>10% of the land area, based on data 
aggregated from FAO [http://www.fao.org/nr/water/aquastat/irrigationmap/index.stm]), 
suggesting difficulty in relating observations-based FPAR values to local precipitation 
amounts. 

Figure 8. (a) Same as Figure 2f (average model-simulated FPAR in peak FPAR season as a 

function of the mean precipitation and the standard deviation of annual precipitation), but 
for the case in which the entire globe is forced to be covered by grassland, (b) Same as 
(a), but for the “all tree” case, (c) Same as Figure 6c (product of the interannual variance 
of model-simulated FPAR averaged over the maximum FPAR season and the square of 
the correlation between NDVI and annual precipitation, multiplied by 500), but for case 
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where the entire globe is foreed to be eovered by grassland, (d) Same as (e), but for the 
“all tree” case. 

Figure 9. Global distribution of annual gross primary productivity (GPP, in grams C / m -day) 
for: (a) the control simulation; (b) the simulation in which all land is covered with a 
grassland vegetation type; and (c) the simulation in which all land is covered with a 
deciduous tree vegetation type. 

Figure 10. (a) Same as Figure 4, but for GPP (rather than NDVI or FPAR) produced in the 
control simulation, (b) Same, but for GPP produced in the “all grass” simulation, (c) 
Same, but for GPP produced in the “all trees” simulation. 

Figure 11. (a) Average GPP (g Carbon / m -day) as a function of the mean precipitation (x- 
coordinate) and the standard deviation of annual precipitation (y-coordinate) in the 
control simulation, with individual land grid cell values averaged over bins. At least 5 
dots must lie within a bin for the binned value to be plotted, (b) Same, but for the 
simulation in which all land is covered with a grassland vegetation type, (c) Same, but 
for the simulation in which all land is covered with a deciduous tree vegetation type. 

Figure 12. Difference in the mean annual GPP produced in the ClimP simulation (the simulation 
using climatological precipitation forcing) and that produced in the control simulation, in 
units of grams carbon/m -day. 

2 4 2 

Figure 13. (a) Variance of annual GPP (in g / m day ) as produced by the control simulation. 

(b) Corr (GPP,P), i.e., the fraction of the GPP variance associated with interannual 
variance in annual precipitation, (c) Same as (b), but for the fraction of the GPP 
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variance associated with interannual variations in annual temperature, (d) Same as (b), 
but for the fraction of the GPP variance associated with variations in annual net radiation. 

Figure 14. Interannual variance of GPP (g /m day ) for the 3-month period centered on the 

month for which the local GPP is climatologically largest. Overlain on the plot are black 
dots showing where the interannual variance of 3 -month FPAR is maximized. 
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Figure 1. a. Schematic of flux computations in the original Catchment LSM. b. Schematic of 
flux computations in the merged model, Catchment-CN. 


747 


748 


36 








749 

750 

751 

752 

753 

754 

755 

756 

757 

758 

759 

760 


Valley bottom 


hillslope 


hilltop 


Dynamic 

hydrological 

zones 

(areas change 
with time) 


ARl 


W, 


AR2 


W, 


Wv2 = weighted average contribution 


from and 


AR3 


W, 


Static 

carbon 

zones 

(areas fixed) 



VI 


V2 


V4 


Figure 2. Schematic of independent vegetation (carbon) treatments in topographically-defined 
vegetation zones. Three static vegetation zones are defined, with independent carbon prognostic 
variables in each. Wi, W2, and W3 are soil moisture states in the three dynamically-varying 
hydrological zones (with time-varying areas ARl, AR2, and AR3); weighted averages of these 
states (e.g, Wv 2 , as shown in the figure) and corresponding weighed temperature states are 
passed down to the carbon physics calculations for the different vegetation zones. Weighted 
averages of vegetation zone quantities (e.g., canopy conductance) are similarly passed back to 
the hydrological zones. 
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763 Figure 3. a. Distribution of average NDVI in peak NDVI season (the month for whieh the mean 

764 seasonal eycle of NDVI is maximized along with the preceding and following months), from 

765 GIMMS observations. The peak season varies with grid cell; see text for details, b. Same, but 

766 for average FPAR in peak FPAR season, from GIMMS observations, c. Same, but for average 

767 FPAR in peak FPAR season, from model simulation. 
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Figure 4. (a) Square of the spatial eorrelation eoeffieient (over land grid eells) between observed 
NDVI and land surface forcing variables: the logarithm of the mean annual precipitation (P), the 
logarithm of the standard deviation of annual precipitation (op), mean annual air temperature (T), 
and mean annual net radiation (Rnet)- The final two bars show the square of the correlation 
coefficient obtained from the multiple regression of NDVI against, respectively, (i) precipitation 
mean and standard deviation, and (ii) all four quantities, (b) Same, but for observed FPAR. (c) 
Same, but for modeled FPAR. 
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779 Figure 5. (a) Average GIMMS NDVI in peak NDVI season as a function of the mean 

780 precipitation (x-coordinate) and the standard deviation of annual precipitation (y-coordinate). 

781 Each dot represents a single land grid cell, (b) Same as (a), but with the individual values in the 

782 scatter plot averaged over bins. At least 5 dots must lie within a bin for the binned value to be 

783 plotted, (c) and (d): Same as (a) and (b), but for GIMSS FPAR data, (e) and (f): Same as (a) and 

784 (b) , but for model-simulated FPAR data. 
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c. Var(FPAR)* (multiplied by 500): Model 
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786 Figure 6. (a) Product of the interannual variance of GIMMS NDVI data averaged over the 

787 maximum NDVI season and the square of the correlation between NDVI and annual 

788 precipitation, multiplied by 1000. (b) Same, but for GIMMS FPAR data, (c) Same, but for 

789 model-generated FPAR data, and with the scaling factor changed to 500. 
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791 

792 Figure 7. Map of grassland locations, as derived from the distributions used in the Second Phase 

793 of the Global Soil Wetness Project (Dirmeyer et al. 2006). The dots overlain on the plot indicate 

794 regions for which irrigation is extensive (>10% of the land area, based on data aggregated from 

795 FAO [http://www.fao. org/nr/water/aquastat/irrigationmap/index.stm]), suggesting difficulty in 

796 relating observations-based FPAR values to local precipitation amounts. 
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800 Figure 8. (a) Same as Figure 2f (average model-simulated FPAR in peak FPAR season as a 

801 funetion of the mean precipitation and the standard deviation of annual precipitation), but for the 

802 case in which the entire globe is forced to be covered by grassland, (b) Same as (a), but for the 

803 “all tree” case, (c) Same as Figure 6c (product of the interannual variance of model-simulated 

804 FPAR averaged over the maximum FPAR season and the square of the correlation between 

805 NDVI and annual precipitation, multiplied by 500), but for case where the entire globe is forced 

806 to be covered by grassland, (d) Same as (c), but for the “all tree” case. 



Var(FPAR)* (multiplied by 500): Model (all trees) 


120W 


60W 


eoE 


120E 


43 


a. Mean GPP (g carbon / m^-day); Model (control) 
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807 

808 Figure 9. Global distribution of annual gross primary productivity (GPP, in grams C / m -day) 

809 for: (a) the control simulation; (b) the simulation in which all land is covered with a grassland 

810 vegetation type; and (c) the simulation in which all land is covered with a deciduous tree 

811 vegetation type. 
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814 Figure 10. (a) Same as Figure 4, but for GPP (rather than NDVI or FPAR) produced in the 

815 control simulation, (b) Same, but for GPP produced in the “all grass” simulation, (c) Same, but 

816 for GPP produced in the “all trees” simulation. 


45 



817 

818 

819 

820 
821 
822 
823 


1.0 


ap 

(mm/day) 


0.1 


0 


GPP: Model (all grass) 



GPP: Model (control) 



.1 1.0 10.0 


GPP: Model (all trees) 


ap 

(mm/day) 



1.0 10.0 

Mean P (mm/day) 


2 

Figure 11. (a) Average GPP (g Carbon / m -day) as a funetion of the mean preeipitation (x- 
eoordinate) and the standard deviation of annual precipitation (y-coordinate) in the control 
simulation, with individual land grid cell values averaged over bins. At least 5 dots must he 
within a bin for the binned value to be plotted, (b) Same, but for the simulation in which all land 
is covered with a grassland vegetation type, (c) Same, but for the simulation in which all land is 
covered with a deciduous tree vegetation type. 
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830 Figure 12. Difference in the mean annual GPP produced in the ClimP simulation (the simulation 

831 using climatological precipitation forcing) and that produced in the control simulation, in units of 

832 grams carbon/m -day. 
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c. Cori^(GPP,T) 
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835 Figure 13. (a) Variance of annual GPP (in g / m day ) as produced by the control simulation. 

836 (b) Corr (GPP,P), i.e., the fraction of the GPP variance associated with interannual variance in 

837 annual precipitation, (c) Same as (b), but for the fraction of the GPP variance associated with 

838 interannual variations in annual temperature, (d) Same as (b), but for the fraction of the GPP 

839 variance associated with variations in annual net radiation. 


840 


841 


48 


842 


Var(GPP): Control 



60N 


SON 


SOS 


60S 


180 


120W 


60W 


60E 


120E 


180 


I 


1.00 

0.81 

0.64 

0.49 

0.36 

0.25 

0.16 

0.09 

0.04 

0.01 

0.00 


843 

2 4 2 

844 Figure 14. Interannual variance of GPP (g /m day ) for the 3 -month period centered on the 

845 month for which the local GPP is climatologically largest. Overlain on the plot are black dots 

846 showing where the interannual variance of 3-month FPAR is maximized. 
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